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Abstract 

We consider the time-dependent Gross-Pitaevskii equation describing the dynamics of rotat¬ 
ing Bose-Einstein condensates and its discretization with the finite element method. We analyze 
a mass conserving Crank-Nicolson-type discretization and prove corresponding a priori error es¬ 
timates with respect to the maximum norm in time and the L 2 - and energy-norm in space. The 
estimates show that we obtain optimal convergence rates under the assumption of additional 
regularity for the solution to the Gross-Pitaevskii equation. We demonstrate the performance of 
the method in numerical experiments. 


1 Introduction 

When a dilute gas of a certain type of Bosons is trapped by a potential and afterwards cooled down 
to extremely low temperatures close to the absolute minimum of 0 Kelvin, a so called Bose-Einstein 
condensate (BEC) is formed [T91I2H 27, 50]. Such a condensate consists of particles that occupy the 
same quantum state. That means that they are no more distinguishable from each other and that 
they behave in their collective like one single ’super-atom’. Recent overviews on the mathematics 
for Bose-Einstein condensates are given in mi ma¬ 
in this work, we focus on the specific case of Bose-Einstein condensates in a rotational frame 
[29 ]. One of the interesting features of a Bose-Einstein condensate is its superfluid behavior. In 
order to distinguish a superfluid from a normal fluid at the quantum level, one needs to verify 
the formation of vortices with a quantized circulation (cf. [2] for an introduction in the context 
of BECs). In experimental setups the formation of such vortices may be triggered by rotating the 
condensate. This can be achieved by using a stirring potential which is generated by imposing 
laser beams on the magnetic trap (cf. [H 051 0U .55:, SU E]). If the rotational speed is sufficiently 
large, the vortices can be detected (cf. [1]). In particular, the equilibrium velocity of the BEC 
can no longer be identified with a solid body rotation and it can be observed that the rotational 
symmetry breaks (cf. |53j for an analytical proof). The number of vortices strongly depends on 
the rotation frequency. However, if the rotational speed is too low no vortices arise and if the 
rotational speed is too large (relative to the strength of the trapping potential) the BEC can be 
destroyed by centrifugal forces. Analytical results concerning the formation, or lack, of vortices, 
their stability, types and structures depending on the rotational speeds and trapping potentials is 
found in [31 [1211231 El El [53]. Detailed numerical investigations are given in Hang. 

The formation and the dynamics of BECs are typically modeled by the Gross-Pitaevskii equation 
(GPE) which is a Schrodinger equation with an additional nonlinear term that accounts for particle- 
particle interactions [3311421149j . To account for a rotating BEC, it is common to extend this model 
by an angular momentum term. Let D C M rf , d = 2, 3, be a bounded convex Lipschitz domain and 
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[0, T] Cla time interval. We consider the dimensionless time-dependent Gross-Pitaevskii equation. 
For the case d = 3 we seek the complex-valued wave function m:Dx[ 0,T]->C that describes the 
quantum state of the condensate. It is the solution with initial state u(-, 0) = uq to the nonlinear 
Schrodinger equation 

i dtu = —-A u + V u + if2 • (x x V) u + /3\u\ 2 u in V, (1) 

u = 0 on dV, 

where we denote x = ( x,y,z ) E M 3 . Here, V characterizes the magnetic trapping potential that 
confines the system (by adjusting V to some trap frequencies) and the nonlinear term (3\u\ 2 u de¬ 
scribes the species of the bosons and how they interact with each other. In particular, /3 depends 
on the number of bosons, their individual mass and their scattering length. We assume that /3 
is strictly positive (which means that we have a repulsive interaction between the particles). The 
term iSI • (x x V) m characterizes the angular rotation of the condensate, where fl E M 3 defines the 
angular velocity. As usual, the operator L = (C x , C y , C z ) := —i(x x V) = x x P describes the 
angular momentum, with P = —iV denoting the momentum operator. 

In the following, we assume that the rotation is around the z-axis, which leads to the simplifica¬ 
tion ifi • (x x V) = —QC Z , where C z = —i ( xd y — yd x ) is the ^-component of the angular momentum. 
With this simplification the weak formulation of problem (fTl) (respectively its dimension reduced 
version in 2d) reads: find u E C°([0, T), Hq(T>)) and dtu E Cq[0, T), H~ 1 (V)) such that u(-, 0) = uq 
and 

i{dtu(-,t), (f>) L 2 (p) = ^{Vu{-,t),V(f )) La(X) ) (2) 

+{Vu(-,t),(f>)L 2 (p) - Sl{£ z u(-,t),4>) L 2 (d) +/3{\u(-,t)\ 2 u(-,t),(j)} L 2( D) 

for all <fi E Hq(T>) and almost every t £ (0, T). Here, (•, -)l 2 (v) denotes the standard L 2 -scalar 
product for complex valued functions, i.e. (v, to)z, 2 (i>) = fp v(x)w(x) dx for v, w E L 2 (V). 

A recent existence and uniqueness result concerning the solution of Q can be found in [8j for the 
case of the three dimensional Cauchy problem, i.e. for the case T> = M 3 ^(see also [36j for an earlier 
work). A general comprehensive overview on existence and uniqueness of nonlinear Schrodinger 
equations can be found in the book by Cazenave [2?] . 

The literature on the numerical treatment of ([2]) is rather limited for the case fl ^ 0. Very 
efficient methods that exploit Fourier expansions were proposed in pa na nm : in DEI a time¬ 
splitting method is proposed that is based on the scaled generalized-Laguerre, Fourier and Hermite 
functions, whereas in m it is suggested to discretize ([2]) in rotating Lagrangian coordinates. A finite 
difference discretization is discussed in |12j. A comparative overview on different time-discretization 
is given in [6]. Concerning the numerical treatment of the eigenvalue problem associated with <©> 
we refer to mm- 

Even though spectral and pseudo-spectral methods (such as the explicit methods proposed in 
mmm) are typically computationally cheaper than a pure finite element based approach as 
proposed in this paper, they generally require a high smoothness of the magnetic potential to 
work. Non-smooth potentials can for instance arise in the context of investigating Josephson effects 
(cf. [59j 60]) or experiments involving very rough disorder potentials (cf. @8]). In corresponding 
numerical simulations, the usage of finite elements seems to be unavoidable for an efficient method. 
Another advantage of finite elements is that they can be easily combined with mesh adaptivity, as 
it is helpful to resolve localized vortices. 

There are some results concerning the convergence of numerical methods with respect to the 
space discretization. Concerning PI finite elements and for the particular case that V = 0, Q = 0 
and that the spatial mesh is quasi-uniform, a priori error estimates can be found in |5j, (38l [39 , 52l 
EH ESI EE]- We describe those results in chronological order. 
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The first results were obtained by Sanz-Serna |52] who considered a modified Crank-Nicolson 
scheme that conserves the mass and the energy. For the case cL = 1 and a periodic boundary 
condition, optimal L 2 -error estimates were derived with a quadratic order convergence in time. A 
necessary condition for the analysis was that the time step size r can be bounded by the mesh size 
h, i.e. the time step size is constrained by r = 0(h). 

In [5j, Akrivis et al. generalize the a priori L 2 -error estimates of [52] to d = 1,2,3 and to the 
case of a homogenous Dirichlet boundary condition. Furthermore they could relax the constraint 
for the time step size to the condition t = 0(h d / 4 ). Beside the modified Crank-Nicolson scheme, 
the authors also study a one-stage Gauss-Legendre implicit Runge-Kutta scheme (IRK) that we 
will also consider in this paper. The IRK is still mass-conservative, but does no longer conserve the 
energy. However, as we will see numerically, the energy deviation is marginal. Again, the condition 
r = 0(h d / 4 ) is required. Furthermore, the authors propose and analyze a Newton-scheme for solving 
the nonlinear problems that arise in each time step. 

In [57], Tourigny investigates the case of optimal L°°- and iH 1 -error estimates. He analyzes the 
same implicit Runge-Kutta scheme as considered in [5] and recovers the constraint r = 0(h d / 4 ). 
Furthermore, he investigates a classical Backward-Euler discretization for which the more severe 
constraint r = 0(h d / 2 ) is required. However, as we will see in our analysis below, both constraints 
are not optimal. For instance for the Backward-Euler scheme, for any s > 1, we can improve it to 
t = C>(| ln/i| -s / 2 ) for d = 2 and to r = 0(h s / 2 ) for d = 3 (see Theorem 


3.5 


below). 


Concerning higher order schemes (without conservation properties), a space-time finite element 
method was proposed and analyzed in [38, (39| for the case d = 2 and for graded meshes. Here, 
[38] is devoted to the case of a Discontinuous Galerkin time discretization and [39] is devoted to 
a Continuous Galerkin time discretization. Optimal error estimates in L 2 and H 1 are derived. 
Here the constraint (for d = 2) reads t p = 0(\ ln/i|^ s ) for some s > 1 and where p denotes the 
polynomial degree used for the time discretization. Hence, it excludes lowest order schemes such as 
the Backward-Euler scheme for which p = 0. 

In ]dl| , Zouraris considers a mass conservative linearly implicit two-step finite element method. 
Zouraris proves optimal order L 2 - and i/ 1 -error estimates under the mild time step conditions 
t = 0( | ln/i| -1 / 3 ) for d = 2 and r = (^(/i 1 / 3 ) for d = 3. 

In a recent work [58] . Wang studies a new type of a linearized Crank-Nicolson discretization 
which is mass but not energy conservative. Again, optimal order L 2 -error estimates are derived 
however with the breakthrough that no constraints for the coupling between time step size and 
mesh size are required. The condition of quasi-uniformity is still necessary. 

Concerning the convergence of space discretizations for the nonlinear GPE eigenvalue problem 
(again for H = 0) we refer to [21] for optimal convergence rates in Fourier and finite elements spaces 
and to m for a two level discretization technique based on suitable orthogonal decompositions. 
Regarding the Gross-Pitaevskii equation with rotation term (i.e. 11 / 0), we are only aware of the 
work by Bao and Cai [12] where optimal error estimates for the finite difference method are proved. 
So far, there seem to be no results concerning finite element approximations. 

In this work we present an error analysis for a Crank-Nicolson-type finite element approximation 
of the time-dependent GPE with rotation. More precisely, we analyze the one-stage Gauss-Legendre 
implicit Runge-Kutta scheme earlier considered by Akrivis et al. |5j and Tourigny m ■ We gen¬ 
eralize these work with respect to two points: we consider the equation with potential and with 
an angular momentum rotation, and for arbitrary s > 1 we show that the time step constrained 
t = 0{h d / 4 ) can be relaxed to r = 0{ | lnfi| -s / 4 ) for d = 2 and to r = 0(h s / 4 ) for d = 3. We do not 
consider Fourier approaches here (even though they can be computationally more efficient in many 
applications), since they require smoothness of the trapping potential, whereas the strength of finite 
element approaches lies in the fact that they do not require such smoothness and that it can be 
easily combined with adaptive mesh refinement strategies. This might be necessary in experiments 
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involving disorder potentials. 

Outline. In Section [2] we establish our model problem and state the basic preliminaries. The 
main results are presented in Section [3j where we state a Crank-Nicolson-type time- and PI Finite 
Element space discretization of the Gross-Pitaevskii equation. Furthermore, corresponding a priori 
error estimates in the L°°(L 2 )-norm and in the L°°(I4 1 )-norm are given. The proof of these estimates 
takes place in several steps. First we introduce a general framework and some auxiliary results by 
investigating the fully continuous problem in weak formulation. This is done in Section [4j In Section 
[5]we show well-posedness of the numerical scheme presented in Section|3j Furthermore, we introduce 
a regularized discrete auxiliary problem which will turn out to produce the same solutions as the 
considered Crank-Nicolson-type Finite Elemente scheme (under suitable assumptions). Finally, in 
Section [6] we derive an error identity and estimate the arising terms. At the end of this section, all 
results are combined to finish the proof of the main theorem. We conclude the paper with numerical 
experiments in Section [7j 


2 Model problem and preliminaries 

Let d = 2,3 denote the space dimension. In order to keep our analysis as general as possible, we 
subsequently consider a slightly generalized Gross-Pitaevskii model. Before stating the problem and 
a corresponding set of assumptions, we introduce our basic notation. 

By x we denote the complex conjugate of a complex number x E C, by x • y we denote the 
Euclidean scalar product between x, y E C d (i.e. x y := Yli=i x iVi) an d by |x| := -y/x • x we denote 
the corresponding norm. The real part of a complex number is denoted by 5ft and by 9 its imaginary 
part. We furthermore use the standard notation for the Sobolev spaces W k,p (T>) (for 0 < k < oo 
and 1 < p < oo) equipped with the norm 

\\v\\ w ,. r{v) := ( ( E WS* l|S “' ,|l W)) 1/ ' ’ for 1 £ P < °°’ 

^max| a |<fc \\d a v\\ L oo(p), for p = oo. 

For p = 2 we write as usual H k (T >) := W k,2 (T>). The semi-norms on H k (T >) are denoted by 

M H*(V) ■= ( X] 

\l«l =k / 

We consider the following model problem. 

Definition 2.1 (Model problem). We consider the (smooth) linear differential operator L : Hq(T>) —> 
that is associated with the following bilinear form, 


{L(v),w) H -np^ H i(py.= / A(x)Vu(x) • Vio(x) + ift(x) • Vv(x)w(x) + c(x)v(x)w(x) dx. (3) 

Jv 

With the above definition we seek u € L°°([0, T), Hq (V)) and dtu € L°°([0, T), H~ l (V)) such that 
u(-, 0) = uq and 


i(d t u{-, t), w) L i (v) = (L(u(-,t)),w) H -i(p )>H i(p) + ((«(■) + P\u(-, t)\ 2 )u(-, t),w) L 2 (v) (4) 

for all w E H ( \ (T>) and almost every t E (0,T). Note that any such solution automatically fulfills 
u E G°([0, T), L 2 (V)) so that u(-, 0) = ilq makes sense. 

Here we make the following assumptions. 
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(Al) The computational domain DcR^ (for d = 2, 3) is a convex bounded polyhedron. 

(A2) The coefficients A, b and c are real valued, smooth and bounded (i.e. L represents the smooth 
linear part of the problem). On the other hand, we assume k E L 00 (T ) , C) H W 1,3 (T>, C); 
P E M > 0 and u 0 E H 2 (V ) D H$(V). 

(A3) The real matrix-valued coefficient A = A(x) is symmetric and there exist positive constants 
7min > 0 and j max > 7 min such that 

7min|£| 2 < A(x)£ ■ £ < 7max|?| 2 for all (£, x) E R d x V. (5) 

By the properties of A there exists a pointwise invertible matrix-valued coefficient A 1 / 2 such 
that A 1 / 2 A 1 / 2 = A. We denote its inverse by A -1 / 2 := (A 1 / 2 ) -1 . 

(A4) The real vector-valued coefficient b = 6 (x) is divergence free, i.e. V • b = 0. 

(A5) It holds 3 ?(k) > 0 and the real coefficient c = c(x) is such that there exist real-valued constants 
Co > 0 and Ci > 1 with 

4c(x) — (2 + Ci)|A -1 / 2 (x)6(x)| 2 > 4Co > 0 for all xgP, 

We note that assumptions (A1)-(A4) are trivially fulfilled for the Gross-Pitevskii equation <©• 
In practice, A is typically just a constant, whereas b describes the angular momentum rotation 
like in (J2]). The term c describes any kind of real-valued non-negative smooth potential such as 
harmonic potentials of the structure c(x) = 'y 2 x 2 + 7 2 y 2 + ^ 2 z 2 with scaled trapping frequencies 
7 au 7 ?/) 7 z £ M. The coefficient k can be used to model arrays of quantum wells for investigating 
Josephson oscillations (see |591 60]) or any other type of rough potential. Furthermore, k can be 
also used to describe imaginary potentials, see for instance the complex double-well potential in [30J 
or applications with phenomenological damping. 

Observe that (A4) implies that the operator L is self-adjoint. Assumption (A5) is an additional 
(often crucial) physical constraint, which says that the rotational speed 0 should be balanced by 
the trapping potential V in the sense that V — 3/2|0 | 2 (x 2 + y 2 ) > Co > 0 on V. The physical 
interpretation is that the trapping potential should be stronger than the arising centrifugal forces. 
Otherwise particles can escape from the trap and the Bose-Einstein condensate is destroyed (hence 
there exist no physically meaningful solutions). As we will see later, the differential operator L is 
elliptic, but degenerates for the case Co = 0 and Ci = 1, which just resembles the instability. 

Remark 2.2. Observe that assumption (A5) allows to balance c and k in a suitable way. For in¬ 
stance, if we only have c > 0 but 49?(k) — (2+Ci)|A -1 / 2 6| 2 > 4Co > 0, we can define «7 ew (x) := k(x) — 
4 -1 (2 + Ci)|A - 1 / 2 (x) 6 (x )| 2 — Co and accordingly c new (x) := c(x) + Co + 4 -1 (2 + Ci)|A _ 1 / 2 (x) 6 (x)| 2 , 
which again suit our assumptions above. Also note that we can hide any imaginary part of c in k 
(which is allowed to be imaginary without constraints). 

Remark 2.3 (Existence and uniqueness). In the case A(x) = 1, 6 (x) = 0, c(x) = 0, k E L°°(X>,M) 
and uq E H ( \ (T>). equation @ admits at least one solution for any time T > 0. The corresponding 
results can be e.g. found in [221 Theorem 3.4.1, Corollary 3.4.2], If d = 2 the solution is also 
unique (cf. [ 221 Corollary 4.3.3 and Remark 3.6.4]). Even though we are not aware of an explicit 
result that guarantees existence of a solution to problem Q under the more general assumptions 
(Al)-(A5), it appears straightforward by exploiting Galerkin’s method and compactness results via 
energy conservation (cf. [22] or [28|, Chapter 7.1 and 7.2]). 
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3 Discretization and main result 


In this section we propose a space-time discretization of problem Q and we state corresponding a 
priori error estimates in L°°(L 2 ) and L°°(H l ). 


3.1 Space discretization 

In the following, we denote by Th a conforming family of partitions of V C R rf that consists of 
simplicial elements and which are shape regular, i.e. there exists an /i-independent shape regularity 
parameter p > 0 such that (for all Th) it holds 

diam (Bk) > pdiam(lL) (6) 

for all I\ E Th, where Bk denotes the largest ball contained in I\. The diameter of an element K E Th 
is denoted by hx', the maximum diameter by h max := maxx<=T h and the minimum diameter 
by h min := minftfeTff ^K- Finally, by h : V M>o we denote the corresponding mesh function 
with h(x) := hx if x E K. For brevity, we subsequently write \\hv\\ H k^ for some v E H k (T>) to 

abbreviate yYh K &T h h 2 K \\v\\ 2 Hk ^ K A . The considered PI Lagrange finite element space Sh C Hq(V) 
is given by 

S h := {v E Hl(V) | VK E T h , v\k is a complex-valued polynomial of total degree < 1}. (7) 

By {Ai,..., Xn h } we denote an ordered (Lagrange) basis of Sh- In particular, we denote by 
Nh =dim(S'/j) the number of degrees of freedom in Sh (which is twice the number of interior nodes 
in Th)- On Sh, we introduce the corresponding L 2 -projection and the Ritz-projection associated 
with L. 

Definition 3.1 (L 2 -projection). The L 2 -projection P L i : Hq(T>) —> Sh is given by 

for v E Hq(T>) : (P L ^(v), w h ) L i(v) = (v, w h ) L 2(p) for all w h e S h - 

Definition 3.2 (Ritz projection). For v E Hq(D) the Ritz-projection Ph(v) £ Sh associated with 
L is given as the unique solution to the problem 


(L(v - Ph(v)), w h ) H -i(p) tH i(p) = 0 for all w h E S h - 


Existence and uniqueness of Ph( v ) follow from Conclusion 4.3 below. 


( 8 ) 


In order to derive the final a priori error estimates, we require further assumptions on the grid 
Th, which will be posed indirectly in the following way exploiting the projections. 


(A6) We assume that the L 2 -projection is IL 1 -stable, i.e. there exists a h-independent constant 2 
such that 

\\ P L^{y)\\ H i{v) < C L 2 IMIhrd) for all v E (V). (9) 

(A7) For p > d, we assume that the Ritz projection given by ([8]) is IT 1,00 -stable for functions in 
W 2 ’^(V), i.e. there exists a /i-independent constant C w i,00 such that 

|| ^Ph{w) 11loo( 25) < C'wl.oo || Vu»|| L co (2 j) (10) 

for all w E Hq{T>) n W 2 ^(T>). Note that since T> is a convex domain, we have the embedding 
W 2 ^{V) w W 1,0 °(V). 
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Both assumptions (A6) and (A7) can be fulfilled by making suitable assumptions on Th■ In 
this paper we directly assume stability of the projections to avoid complicated mesh assumptions. 
Concerning (A6), recent results on the T/Cstability of Pj 2 on adaptively refined grids can be found 
in pirn E]. Concerning (A7), we refer to [2Q; Theorem 8.1.11] where the result is established for 
quasi-uniform meshes. For results on graded meshes we refer to [33, 20 • We note that, the results 
on graded (locally quasi-uniform) meshes are only proved for the Laplacian operator, i.e. L = —A, 
and its generalization to general elliptic operators is still open. However, it seems to be crucial that 
the operator L is sufficiently smooth for (A7) to hold on graded meshes, this is why it might be 
important that k is not included in L (in this context, see also the Holder-estimates for the Green’s 
functions proved in [34] and the necessary regularity assumptions made in [37]). 


3.2 Time discretization, method and main result 

In this paper we assume that the time interval [0,T] is divided into 0 =: to < H < • ■ ■ < Gv : = T. 
Accordingly we define the n’th time interval by I n := t n \, the n’th time step by r n := t n — t n -1 

and step size function t G L°°(0,T) by tu := T n . For simplicity we subsequently only write 
(-, •) := (■, ■)h^ 1 (V),h ip) for the dual pairing on i7 1 (P). We consider the following one-stage Gauss- 
Legendre implicit Runge-Kutta scheme (which is of Crank-Nicolson-type). The scheme is mass 
conservative provided that 9 ?(k) = 0. 

Definition 3.3 (IRK Method for GPE). Let u® := Ih( u o) G Sh be the Lagrange interpolation of 
uq. For n > 1, we seek the approximation u]) G Sh with 

(v>h> v h)tf(p)+T n i{L(ul 2 ),v h ) +r„i((K + /3|«J 2 | 2 )u n h 2 , v h } L 2 {v) = « _1 , v h ) L 2 {v) (11) 


for all Vh G Sh and where u. 


# + < _1 


)/ 2 - 


For alternative time-discretizations based on operator splitting for nonlinear Schrodinger equa¬ 
tions with a cubic nonlinearity we refer to [431 [32] . for the case without rotation, and to [6], for the 
case with rotation. More general approaches are discussed in [35] . 

We note that the IRK scheme given by © is mass conservative if Q(k) = 0. The mass 

conservation, i.e. ||u](||x, 2 (x>) = | I'M/I 11 £ 2 CD) f° r n — 0> is immediately seen by testing with u ^ 2 
in © and taking the real part. Note that the conservation property implies that the scheme is 
unconditionally Testable. 


The following proposition follows from Lemma |5.5| and the proof of Theorem 3.5 below. 


Proposition 3.4 (Existence and uniqueness). If (Al)-(A5) are fulfilled and if h and r n are small 
enough and such that 4(li max + r^) —> 0 for h,r n —> 0, then there exists a solution u^ of ©• If 
Q(k) = 0 and if r n is sufficiently small compared to h (in the sense of Lemma 5.5 below), then the 
solution is also unique. 


The main result of the work is the following a priori error estimate, which we prove in Section 
[6j Recall assumptions (A1)-(A5) from Section [2] and (A6)-(A7) from Section 

Theorem 3.5 (Error estimates for the IRK discretization). Let assumptions (Al)-(A7) be fulfilled, 
let u G IF 2,oo (0,T; H 3 (V)) denote a solution of &, and let h and r n be such that Ihf h ™^ + t 2 ) —> 0 
for h, r n —> 0 and where 


3.1 


In h min I 1 / 2 


for d = 2 
for d = 3. 
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Then, if h and r n are small enough, there exist generic constants C = C(u) that are independent of 
h, r n and T such that for a solution u / of (3.3) and for m E {0,1} it holds 

\\u(-,T) — < C\h 2 m u(-,T)\ H 2 ( V -) + Ce CT (|/i 2 m «o|_ff2(x>) + ll^ 2 m 9tu\\L 2 (o,T,H 2 (v))) 


+Ce J 



1/2 


2 —m 


U 


ll°°(4,fi' 2 (D)) + ll' r2 ^ti' u llioo(j fc) ij2+, 


, || 2 || 2 

l (X>)) + H r “II W 2 ’°°{I k ,H™(T>)) 


We observe that the method yields optimal convergence rates, i.e. it is of quadratic order in 
space and time for the L 2 -error and of linear order in space for the H 1 -error. Details on the arising 
constants in Theorem 13.51 can be found in Lemma 16.61 and Lemma 16.71 below. 


Remark 3.6. It is surprising that the L°°(H 1 )-e stimate in Theorem 3.5 requires the higher regular¬ 
ity dttu(-,t) E H^ifD). A similar observation has already been made by Karakashian and Makridakis 
[38l Remark 4.3] for the simpler equation i dtu = —A u + f3\u\ 2 u. It should be investigated in the 
future if it is possible to weaken this regularity assumption. 

Finally, let us state the corresponding result that can be derived for the Backward-Euler Method. 
This result is rather for comparison, since the Backward-Euler is practically not desirable since it 
lacks both mass and energy conservation. 

Theorem 3.7 (Error estimates for a Backward-Euler discretization). Assume (A1)-(A7), u E 
VF 1,2 (0, T; H^(V)) and h and r n such that lh{h ma .* + T n) ~> 0 for h,r n —>• 0. Let further u ° := 
Zh{u o) E Sh- Then, for all small enough h and r n , there exists E Sh with 

«, v h ) L 2 {D) + T n i ( L{ul),v h ) + T n i ((« + /3|«h| 2 )«h, Vh)^(v) = v h ) L 2 (c) (12) 

for all Vh G Sh and there exist generic constants C = C(u) that are independent of h, T n and T 
such that for m E {0,1} 

\\u(-,T) — Uh\\ H m(p) < C\h 2 m u(-,T)\ H 2 (p) + Ce CT (jh 2 m uo\ H 2 ^ + \\h 2 7n 9iu|| L 2( 0jT ^2(x)))) 


+Ce 


CT 


'2—m 


u \\l 2 (o,t,h 2 (V)) + \\Tdtu(-,t)\\ L 2 ( QiT ' H m+i(p' ) ') + ^2 T k\\h 2 m u(-,t k )\\ h 2 (V) ■ 


k= 1 


The proof of this theorem exploits the same techniques as the one of Theorem 3.5, which is why 
we will not present it here. 


4 Reformulation of the continuous problem 

In this section, we establish some auxiliary results and preliminaries concerning the model problem 
|4]). In particular, we introduce a suitable scalar product on H l (fD) which can be associated with 
the operator L and which is more convenient for the analysis in the following sections. 

If clear from the context, we subsequently leave out the integration variable in our integrals, for 
instance we write v for f D v(x) dx. In order to analyze problem Q properly, we require some 
additional definitions and auxiliary results. 

Definition 4.1. For any subdomain cu C V we define the sesquilinear form (•, -)e(u) by 

{v,w) E ( u y.= J ^A 1 / 2 Vu — i2 -1 A -1 / 2 6u^ • (A 1 / 2 Vrc — i2 _1 A _1 / 2 6m) 

+ [ (c- (1/4)| A~^ 2 b\ 2 )vw 

J (jJ 

for v, w E H 1 (w). Note that c— (1/4)| A~ 1//2 6| 2 is positive by (A5). Accordingly, we define the norm 

II • II e{w) b y IMU(«) : = ^(v,v) E{u) . 
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Lemma 4.2. Let oj C Ll be a subdomain. Under assumptions (A1)-(A5), the sesquilinear form 
(v)-E(uj) is a scalar product on H 1 {uj) and the induced norm |M|_e( w ) equivalent to the standard 
H 1 -norm || ■ ||#i( w \. In particular we have for all v G H l [u) 

ll' y \\e(w) - (! ~ Cl l )\\A l/2 Vv Ili 2 ( w ) + Co|Ml£ 2(w) - 

Proof. Obviously, (-,-)e(uj) is a symmetric sesquilinear form on H ( ] (oj). Hence, it only remains to 
show the existence of constants ce and Ce such that 

< C E \\v\\ 2 h i { cj) for all v G H l {w). 

The upper bound is straightforward using the boundedness of the coefficients. To verify the lower 
bound, we first observe with Youngs inequality for any e > 0 that 

[ | A^Vv-12- 1 A~V 2 bu\ 2 


> 


AX/vVv- f \A~ 1 / 2 b\\v\\A 1 / 2 Vv\ - - f \A~ l/2 b\ 

' J OJ J OJ 


21 |2 
l v\ 


> (1 - e” 1 ) [ AX7v • X7v - 


\A~ 1/2 b\ 2 \v\ 2 . 


Hence 


(v,v) E ^)> (l-e x ) f AVvVv- 2 ^ f \A 1/2 b\ 2 \v\ 2 + 

J OJ ^ J OJ 


C\V\ 


Choosing e = Ci together with (A5) finishes the result (where we assumed Ci > 1). Also observe 
that Ci < 1 leads to degeneracies. □ 

Conclusion 4.3. The differential operator L is uniformly elliptic and continuous on Hq (T>). In 
particular it holds 


(v,w)e(v) = (L(v),w) for all v,w G Hq(U). 


(13) 


Observe that Lemma |4.2| and Conclusion 4.3 imply that the operator L degenerates for Co = 0 
and Ci = 1- 

Proof of Conclusion f.3 . Let v,w G Hq(V), we observe that 

J ^A 1 / 2 Vu — i2 ~ 1 A~ ll,2 bv^ ■ (A 1 / 2 Vrc — i2 ~ 1 A~ 1 / 2 bw) 


AX7v ■ X7w — i2 


-l 


lv 


f vb-X/w — f ro6-Vv)+^- f 

b Jv J 4 J v 


| A 1/2 b\ 2 vw 


AX/v ■ X/w + f wib ■ X7v + \ 

Jv 4 


lv 


where we used that V -6 = 0. Hence we have 


lv 


\A-^ 2 b\ 2 vw, 


,| 2 „ — 


(u, w)e(v) = / AX7v-X7w + / wib ■ X7v + 

Jv Jv 


L 


cvw. 


(14) 


Assumption (A4) finishes the proof of (13). The continuity and ellipticity of L hence follow using 
Lemma 14.21 □ 
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Remark 4.4. Let ui C T> be a subdomain and v,w G i/ 1 (w) arbitrary. Under assumptions (Al)- 
(A5) we see that there exists a constant C (only depending on A, b and c) such that 


AVv ■ Vu> + bS7v ■ w + cvw 


< c , IMIhi (w) II«;||hi ( w ) - 


Using the norm equivalence of Lemma 4.2 we hence also have 


AVu • \7w + bS7v ■ w + cvw 


< C E \\v\\ E (u:)\\w\\e{u)i 


(15) 


with C E = C E (A,b, c). However, note that we do not have {v,w)e(w) = (L(u), w}h-i( lo ),h 1 (u) f° r 
arbitrary v,w G iL 1 (w). 


5 Existence and uniqueness of discrete solutions 

In this section we consider the existence and uniqueness of discrete solutions. For that, we require 
the following well known result which can be found e.g. in the book by Thomee m Lemma 6.4]. 
It can be easily proved using Sobolev embeddings with a inverse inequality. 

Lemma 5.1. Let V C be a convex domain. Then there exists some constant such that for 
all v h G S h 


\ v h\\L°°(V) < CMVVbWw), 


where 


4 : = 


| In h min I 1 / 2 
|4nin| _1/2 


for d = 2 
for d = 3. 


We treat the existence of discrete solutions v% of 0 together with the solutions of some 
regularized auxiliary problem. This auxiliary problem is essential for the analysis of For this 
purpose, we recall a lemma that was basically proved in 

Lemma 5.2. Let M G M be given by 


M \\u ||vy 1 .°°(j„,wu.°°( 2 ))) + CVi,oo(diam(P) + 1) |Vw 

(16) 

where C w i,<x> is the constant from (A 7). Then, there exists a 
cm > 0 such that: 

function fM '■ C — > C and a constant 

/m4) = \z\ 2 z, 

if \z\ < M, 

(17) 

(fM(z),z ) G M>0, 

for all zG C, 

(18) 

\f M (z)\ < 2M 2 \z\, 

for all zGC, 

(19) 

| fM(z) - fM(w) | < 10M 2 \Z - w\, 

for all z,w G C, 

(20) 

II fM(z) - fuiw)\\ E(p) < c M \\z - w\\ E( p) 


(21) 

for all z,w G Hq(V) with 

w Ww 1 ’ °°(T>) — Af. 


The above lemma is a slightly generalized version of |38l Lemma 4.1] in the sense that 
more precise about the constants in (19) and (20), condition (18) is new and condition 

we are 

(21) is 


formulated with a different norm. The latter two points are obvious, therefore we only prove (19) 
and (20). 
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Proof. Let us define 0 := M 2 , g(s) := 39 4 s 5 — 79 3 s 4 + 49 2 s 3 + s and the curve 7 : M —> M by 


70 ) 


s for s < 9 

< g(s — 9) + 9 for s E [6 , 29] 

29 for s > 29. 

— 


It can be verified that 7 E C 2 (M) and we can hence define /m(~) := l(\z\ 2 )z f° r £ £ C. In order to 
verify the properties of /m it is sufficient to check the behavior of g on [0, 9\. It holds 

g'(s) = (150"V + 2 9~ 3 s + 0“ 2 )(s - 9) 2 


which is obviously strictly positive on [0,0). Hence g is monotonically increasing and so is 7 . 
Furthermore, we have 


g"(s) = 600" 4 s(0 - s)(20/5 - s), 


which implies that g' has a maximum in 20/5 with </(20/5) < 2 
Combining these properties of g allows us to derive (19) and d20| 


We observe |</'(s)| < 6O0~ 4 . 
with the constants as given in 
the lemma. Property (18) is obvious since 7 is monotonically increasing and hence non-negative on 

Lemma 4.1] with the //’-seminorm. but follows directly by 

□ 


[0,oo). Condition (21) is stated in 
the norm equivalence that we showed earlier. 


Using the previously introduced function Jm, we can now state the regularized problem. As we 
will see later, the solution to the regularized problem is a solution of the discrete problem © for 
sufficiently small time steps. 


Definition 5.3 (Discrete auxiliary problem). Let /m denote a function with the properties depicted 

E Sh with u° h being the initial value used for problem 
E Sh denote the solution of 


5.2 


in Lemma 

0 - For n > 1 we let U n 


Furthermore we let U° = u° h 


(U U , Vh)L 2 (V) + T n i 


! )) v h) + (kU 71 2 +f3fM(U n 2 ),Vh) L 2 (V)) — {U n 1 ,Vh)L 2 {V) (22) 


for all Vh E Sh and where we defined U n 2 := (U n + U n 1 )/2. 


In order to show existence of the solutions of problem © and \22\ we require the following 
lemma, which is a well-known conclusion from Brouwers fixed point theorem. 


Lemma 5.4. Let N E N and let Bi(0) := {ct E C^| |a| < 1} denote the closed unit disk in C N . 
Then every continuous function g : C N C N with 3?(g(a),a) > 0 for all a E dBi(0) has a zero 
in B i(0), i.e. a point a 0 E -Bi(O) with g(oto) = 0. 

If there exits no e*o E -Bi(O) with g(ao) = 0, then g(a) := — g(a)/\g(a)\ (interpreted as a 
function g : M 2Ar —y M 2Af ) has a fixed point a* E -Bi(O) by Brouwers fixed point theorem. Hence 
1 = \g( a *)\ 2 = (g(a*),a*) = -(g( a *), a *)/\g(<x*)\ = g(a*),at*)/\g(aL*)\ < 0 , which is a 

contradiction. 


Lemma 5.5. For every n > 1 there exists a solution U n E Sh of problem (22). If the time step 
size is such that r n <2 (||Cy(Av)||z.°o( 2 ?) + /310M 2 ) then the solution U n E Sh is also unique. Recall 
that M is the constant appearing in Lemma [5J| Furthermore, if7s(n) = 0 and if r n and h are such 
that r n h~f n —> 0 for r n ,h 0, then the solution u E Sh of problem (11) is unique for sufficiently 


small r n as well. 
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Proof. We start with the existence result for the solution JJ n of problem ( |22[ ). First, recall that 
Nh =dim(jS'/ l ) and that A m denotes the m’th Lagrange basis function. We want to apply Lemma 


N h 

9 z( a ) := — T n i "y ^ a rn (\ m , ^z)l 2 (t>) 

m= 1 

^ Nh i i N h 

+ 2 (f'(Am)i A^) + ((kM + fifM){-jU n 1 + - X()l 2 (V) + -Pi, 

ra=l m=l 


5.4 and define 5 : —>• C Ar, ‘ for c* G C Ar ' 1 by 


where F G C A A is defined by 

Ft ■■= liHU 11 - 1 ), Xi) + (T- l \U n -\\z) L 2 [v) . 


To show the existence of some cxq with g(ao) = 0, it is sufficient (by scaling arguments) to show 
that there exists some K G M>o so that 9?(< 7 ( 0 :), ol) > 0 for all a G C Nh with | a\ = K. For brevity, 
let us denote a := J2m h = 1 Q mA m . Since 9? (r~ l i(a, ch)l 2 (v)) = 0; {a,a)E(v) = (-£(«)> a ) and by 
construction of /m 

/ a + jjn-l , jjn-l \ 

» («rid + Mm)( 2 ), 2 - °’ 

we obtain 

1 11 

sR(g(a),a) > -||a||| (I)) + K(F, cc) - ((/dd +/3/ M )(-U n_1 + 2 X] a ™ A ™)’ un ^)L 2 (v) 

m =1 

\ ^||_||2 ||_—lrm—111 11 _ 11 11 _, 11 \\Tjn—l\\ 

^ 2\\ a \\ E {V)~ W T n U \\L 2 (V)\\ a \\L 2 (V) ~ 2\\ a \\E(V)\\ U II E{V) 

~ 7) (IMIl°°(z>) + (32M 2 ) || U“ 1 ||l2(d) (IMIl 2 (£>) + W 1 ||l 2 (-d)) 

> II«IIb(D) {Ci\\a\\ E( D) - C 2 ) - C 3 , 


where we used the Poincare-Friedrichs inequality in the last step and where C\, C 2 and C 3 are 
appropriate a-independent positive constants. Consequently, for all a with ||cr||^( 7 ?) > C 2 /C 1 + 
\J (C 2 /C 1) 2 + C 3 /C 1 we have 8 t(g(a),a) > 0 and hence, by norm equivalence in finite dimensional 
spaces, there exists a sufficiently large K such that $l(g(a), a) > 0 for all a. with |cc| = K. This 
gives us existence of a discrete solution of ( 22 ). 

For uniqueness in (22) we use an L 2 -contraction argument. Let us compare two solution U ( 
and UJ 1 ^ of problem (22). Using the equation and testing with Ufa — Ufa we get 


(i) 


( 2 ) 

1 ^( 1 ) — ufa\\l 2(v) 


2 1 


i [(L(Ufa - Ufa), Ufa - Ufa) + fa{Ufa - Ufa), Ufa - Ufa ) L2{v) ) 

uj l u + u n ~ 1 


umtu Ufa + U "- 1 Ufa + U n - X Ufa + U "- 1 

-2r n i^(/ M (—- 7 ,-) - /m(— -x-),-x-x- )l 2 (d) 


= -(TnO^XC/fo - t/f 2) ), t/fo - Ufa) L2(D) 


( u n 

(/m( — 




n—1 


U™+U n ~ l U™+U n ~ l Ujl. +u n ~ 1 

) JM\ ry j 5 ry rs / L 2 (T>) 


^ xF> (||9f(«)|| L cc (2?) +/310M 2 ) ||Ufa - Ufa\\ 2 L2{vy 


2 


2 
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Since we assumed that r n (||Qf(/c)]|x,<»(x?) + /310Af 2 ) < 2 we conclude || \\l 2 (d) = 0 and have 
hence uniqueness. 

For two solutions and ^ of the original IRK scheme ( p~I| ) , we use the additional as¬ 

sumption 3 ! (k) = 0 to conclude with the mass conservation that 


\Wh,{l)\\L 2 (V) ~ \\ u h,( 2 )\\L 2 (V) ~ I \ u h II L 2 (V) ~ II w AiIIl 2 (X>) — llA('Uo) ||l 2 (X>) < C\\u 0 \\h 2 (V) > 

where we used the stability estimate \\Th(v 0||l 2 (23) < C1MI-f/ 2 (x>) for v G H 2 {V) for the Lagrange 
interpolation operator Z/,. With this, we can proceed as before to obtain 


l“fc,(l) ~ U h,(2)\\L^V) 


< 

* T i p L 


\(i) 

U h,(l) 
2 


+ K 1 


U 


h,( 2) 


+ < 1 


I n, n ~ D „, n \ 

u h,{ 2) + u h b u /i,(l) _ u h,{2))L 2 (V) 


+<~ 1 


+ 


1 l n 

u h,( 2) 
2 


+ < _1 


■u 


Ai,(l) U h,{2) 


— Ai/?^||«o||tf2(£>) 11^(1) ' u h,(2)lll,° 0 (X>) 


With the inverse estimate \\u^ ^ — u ^ , 2 ^ ||i,oo(X)) < Ch mi ^ 2 ||u^ ^ ^ ||z, 2 (x?) we conclude that for 

an appropriate positive constant 

ll u h,(l) ^ft,(2) Hl 2 (D) — 11 u 0 11 H 2 (V)) r n^milll tt ft,(l) — ^,(2) II i 2 (23) 

and hence uJW, = , 2 , for sufficiently small r n . □ 


6 A priori error estimates 

In the following we assume that u denotes a solution of © with sufficient regularity. In this section 
we derive an a priori error estimate for the discrete solutions. However, instead of taking © as our 
reference problem we follow the ideas of [35] and take the auxiliary problem (22) as our reference. 
In this context, note that by the definitions of u and Jm we have 


(u(-,t n ),v h ) L 2(D)+i/ (L(u),v h )+i {ku +PfM{u),v h ) L 2( V -) = {u{-,t n -i),v h ) L 2(p) ( 23 ) 

J I n J I n 

for all Vh G S^. Since u is continuous in time we can define u n := u(-,t n ). 

For simplicity (and slightly abusing the notation), we write for v G H 2 {T >) D Hq(T>) 


Lv := —V • (AVu) + b ■ Vu + cu, 


so that 


(v,w) E( p) = {Lv,w) l2(v) for all v,w G Hq(T>) D H 2 (V). 


(24) 


In order to derive the a priori error estimates, we first derive an error identity and then estimate 
the various terms in the identity. 

Before starting, recall Definition |3.1[ i.e. the definition of the Ritz-projection associated with 
L. Note that we do not include the term (act;, 4>h) in the Ritz-projection since we want I to be a 


smooth and self-adjoint operator. Since ac can be imaginary, equation (13) would not be valid any 
longer. 

Finally, we also recall a standard result (which follows from the best approximation property of 
Ph with respect to the IL 1 -norm and an Aubin-Nitsche duality argument). 
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Lemma 6.1. Assume (A1)-(A5). There exist generic positive constants C\ and C 2 such that 

lb - p h(v)\\ L 2 (p) < Ci\h 2 v\ H 2 (p) and \\v - Ph(v)\\ E (v) < C 2 \hv\ H 2 ( p ) (25) 

for alive H 2 (V)nH^(V). 

In the first step, we establish an error identity. 

Lemma 6.2 (Error identity). We introduce the abbreviation f(y ) := nv + /3/m( n). For n € N, 
n > 1 we define the error splitting by 


e n h ■.= (Uf - P h {u n )) + (P h (u n ) - u n ) 


and the error contributions by 


d 1} : = Ph{d t u{-fi)) - dtu(-,t) dt, Ci 2) := r n \ (j( 


Phju^ + Ph^- 1 ) _ ~ U n + U n ~ 1 
2 ' ^ 2 


4 3) : = 1 J /(«(■,*)) - /( 


i^(« n ) + ^(« n_1 ) 


) (it, ^ 4) := i J u(-,t ) - 


n n + a ™” 1 


Wit/i these notations the following L 2 -norm identity holds for 

K\\h(v) 

= \\K~ l Whip) + 3* (<d 1} + + $\EZ + + <£(ei 4) ),K + K" 1 )) • 

and the following energy-nom identity 

jKlllio = IK _1 lll(D) + + K~h PlAH i 1 ’ + d 2> + 6 3> » 

+ Si (i(L(Ei + Ep), J' P L 2(Lu(;t)) - MtPPIaAP ‘> dt)) 


Proof. Recalling the definition of U n we have for all € Sh 


Tjn 1 Jjn-Y Tjn _i_ TTn -1 

rn Trn—1 \ _ , ^ ... \ _ +U 


(U n -ljn-\ Vh ) L2{v) + Tn i( L (- 
Subtracting the term 


-)> v fc) = -T-ni(/(- 


)dh)i 2 (C)- 


<^(« n ) - ^(^Wl^) +r Ti ,i(L( PfeK) W 
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on both sides of (29) gives us 


1 


(Eh ~ E n h -\v h ) L . {v) + -r n i (L{E n h + E n h ~ v \v h ) 

= (Ph(u n ~ 1 ),v h ) L 2 (p) - (Ph(u n ),v h ) L i {v) - Tn i (L( 


P h (u n ) + P h (u n - 1 ) 


),V h ) 


n I rrn-1 


~Tni(f( 


U n + U 


-),Vk)l 2 (V) 


= (Ph(u n *), v h } L 2 {v) - (P h (u n ),v h ) L 2 {D) - r n i (L( 


u n + u n 1 


).Wfc> 


+ r„i ( (/( + 


1231) 


=* (Ph(u n *) - u n \ V h ) L Z(V) - (Ph(u n ) - U n , v h } L 2 {V) 


(30) 


u n + u n 1 




+i( / L(u)(-,t) dt — r n L{ 

'I In 

>.,[ h P h (u n ) + P h {u n ~ 1 )^ x 

+!( / f(u(-,t)) dt - T n f{ ---), V h ) L 2( V) 

'I In 

. (, hlPh {u n ) + P h {u n ~ 1 )^ h U n + U n -\ x 

+Tn 1 (/(- 9 -) - /(- 9 -), V h ) L 2 {D) 


Testing with Vh = E% + E^ 1 and only using the real part of the equation gives us 

\\mh{v) 


= i(/( 


'Phju^ + Phju"- 1 ) + ^ | ^_ 1n 


2 , -v 2 -),EZ + E^) L 2 (v) )+\\E" h -% 2{v) 

(( P h (u "- 1 ) - + iOiW - » (W«") - + E^- 1 )^) 

+ » {i(j L(U)(; t ) dt - T n U pPP ‘ ), £J + Ej" 1 )! 

+»(i <1 /(«(■■ 0) <« - Tj{ Ph{u '" > + 2 P '‘ ( “"~ 1) ), EH + iO^(P)) • 

The simplification 

5R (W« n_1 ) - + K _ 1 >^ ( d)) - ^ «^K) - « n , ££ + K _ 1 W)) 

= Ph(dtu(-,t)) — dtu(-,t) dt, Eft + E^~^) L 2(p)^ 

finishes the proof of the L 2 -norm identity. 

To derive the energy-norm identity we use the L 2 -Riesz representer GJJ G of the error func¬ 
tional (L(E%), •). The Riesz representer G G 5/,, is characterized by the equation 


(vfc, = (L(E% + for all G 5 h . 


(31) 




















16 


Testing with Vh = G in (30) and using (L(E% + E'f' ), Gff) = we obtain 

(L(K + K~').K - K~‘) + T„i \\\GIWI HV) 

= (Pftfn"- 1 ) - GDmv) - <«,(«") - G;)y (B) 

,n j „,n— 1 


+i f (L(u)(;t)-L( Ua + ™'‘ ),G n h )dt 

J In 


+i( [ f(u(-,t))dt-r n f( 
Jin 


P h {u n ) + PhiiE- 1 ) 


•> Gh)L 2 (V) 


+ r ni 


2 ' J y 2 
7 <n—l\ zpn Tpn— 1\ _ n rpn\\2 


Using that X(L(E% + i^’ 1 ), - E^ 1 ) = ||^||| (c) - \\E^% {V) and that 

M(v,G%) l2{v) = ^(LiEl + Et- 1 )^^)) for all v € H^(V) 


and taking the real part of equation (32) yields 


!KII|(X>) = ll^r 1 lll (C ) + *t(L(EZ + e^),p l2 (^ + tf) + &) 

+ s (,( i «+ e ;-‘),/ p L ,(Lu(;t))- F ^ (Lvn) + F ^ Lu "" > 

and finishes the proof. 

The next lemma is central for estimating the /-terms in the error identities. 


dt) 


Lemma 6.3. Recall the constant M from (16) and let f(z ) := \z\ 2 z. It holds (a.e. inT>) 

r tn 

/ f(u(-,t))dt-T n f 

'I in — 1 


u n + u"- 1 


<r„Af [ 'J\v n -u n ~ l I 2 


- n n 1 | 2 + T 2 M||u|| W 2 ,oo(/ n )^ 


and 


V 


f(u(-,t )) dt - T n f 


u n + u n 1 ' 


f(u(-,t)) dt - T n f 


’in — 1 


u n + u n 1 ' 


f /(*,*» * - r . ^) + /("- 1) ) + Tn ( /<“■) + / ( " n ~ 1) - / 


i / 1 + it n 1 


(32) 


□ 


(33) 


(34) 


< 4r n M flu" - u "- 1 ! 2 + | V« n - V^" 1 ! 2 ) + r >/ 2 (||«|| w 2 .ao (Jn) + ||V«|| wa ,ao (Jn) ) . 
Proof. We decompose the error under considerations into 


(35) 


With /(it) = \u\ 2 u, the first term in (35) can be estimate using the trapezoidal-rule to obtain 


it it — r, 


\u n \ 2 u n + |it n " 1 | 2 it n - 1 


n 


in— 1 




< T n M \\u\\ W 2,oo( In y 


+ \\ d ttu\\L°°(I n )\\'u\\L°°(I„) 


(36) 
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For the second term in (35), let ( n : [0,1] —> [u n ~ 1 ,u n ] denote the complex valued (linear) curve 
given by Cn(s) := (1 — s)ri n_1 + su n for s € [ 0 , 1 ]. We have Cn( z ) = u n — u n ~ l (and = 0 ). With 
that, we get with the trapezoidal-rule and the midpoint rule that 


f(u n ) + f(u n x ) 


< 


-/ 


(/°Cn)(0) + (/ ° Cn)(1) 


u n + U 


n —1 


[ (/°Cn)(s) 

JO 


ds 


+ 


(f°Cn)(s) ds- f 


Cn (0) + Cn (1) 


< Yo II (/ ° Cn) /, ||L oo (0,l) + ^ll(/ 0 Cn)^ ||l oo (0,1) — g IK/ ° Cn)" ||l°°(0,1) 


12 ' 

. 3 .. 

< -\u — u 

4 


n- 1|2 


IICn||z,°°(0,l) ^ ^ 


^ ° l„,n „n- 1 12m n 

< -|M — U | 


(37) 


can 


Combining the estimates (36) and (37) with (35) finishes the proof of (33). Estimate 
be derived analogously by applying trapezoidal-rule and midpoint rule to the function g(s) := 
2| (1 — s )« n_1 + su n | 2 ((l — s)Vw ’ l_1 + sVtt") + ((1 — s)u n_1 + su n ) 2 ( 1 — s)Vu n ~ 1 + sVu n . □ 

Lemma 6.4 (L 2 -error estimate for EVf). Consider n > 1 and E= U n — Ph{u n ). Let M denote 
the constant in Lemma 5.2 . There exists a constant Cm that only depends on M, V, k, (3, C\y i,oo 

(38) 


and C i (cf. the L 2 -estimate (25)2 such that for all r n < (2 Cm) 1 it holds 


\K\\l 2 { v) - (1 - Cmt \ 11^ X |l L 2 {V) + c M\\h 2 d t u\\ 2 L 2 ( Iri j 1 2 ( T) )) 


+CmTu ^||r 2 L(9 tt u)|| 2 <x> ( Jn)L 2(x))) + \\ T 2 u \\w 2 ’°°(i ri ,L 2 (V)) + \\h 2u \\ l°°(/„,h 2 (x>))) • 

Proof. In the following Cm denotes any constant that depends generically on M, V, n, (3, C w i, 0 
an 
3? 


and C\. We estimate the terms on the right of side of the error identity (26) and start with 
tW 


n ,E% + E t 


n— 1\ 


)L 2 (V) 


. We obtain 


< IlCn, 1 ^ IIl 2 (x>) 


\\E n h +E n h - l W L \v) 

= || f Ph(d t u(-,t)) - d t u(-,t)dt\\ L 2 {D ) < j \\P h (d t u(-,t)) - d t u(-,t)\\ L 2 {D) dt 

J I n J I n 

9 C\ ( \\h 2 d t u(-,t))\\ H 2 (D) dt < CiT^WtfdtuWtf^nz^y 

Jin 


(39) 


Hence 

+ E'f 1 ) L 2(X))) | < 2 Tn (\\EUh iD) + \\E% 1 \\l 2 (V)) + -Cf \\h 2 d t u\\l 2 ( In ' H 2(p)y 

Next we bound the term depending on = r n i ^/( P '^ u )+Pd M -2) _ /( c/n+pn 1 )'j . Recalling 

that = U n — P h {u n ) we obtain 


i3f((ci 2) ,^+ J Er 1 )^p ) )i<^n/( 


P^O + P^rd 1 - 1 ). . c/^ + c/"- 1 .,. . 

-~-) - /(-5-)IIl 2 (d)I 


Ek + E n ~ l 


h^r^h \\l 2 (V) 


9 (|k||Loc (c) + 10M 2 /3)r n (||P^ 


in||2 ill T7' n ~ 1II2 

L 2 (X>) IIHl 2 (d) 


(40) 
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Recall f(z) = nz + / 3f M (z)■ In order to treat ^ 3) = i (f(u) - f (\(P h {u n ) + P h (u n x )) , 1) L 2 (j n) , we 
use that /m(z) = \ z \ 2 z f° r \ z \ < M and the facts that |M|l°°(/„x:d) < M and ||P/ l (u n )||ioo(x>) < 
C' w i,oodiam(R ) )||Vrt n ||/ / oo(-p) < M to conclude that 


l»(«i 3> ,ES + £rW>) 

K + K-'tivm 
+PW [ /(«(■>*))-/( 

■Jin 


.. f , ^ P h {u n ) + ^(« n_1 ) 

< IMIl°°0d)H / “(■)*)-^-^IIlrd) 

J In “ 


iM« n ) + ^(« n_1 ) 


) ^IIl 2 (D)) 


(41) 


where /(z) := To estimate this, we decompose || J 7 f(u(-,t)) — f{ Ph< " u ||z, 2 (x>) 

into 


/(«(•.*)) - / 


+ jj "' 1 


^IIl 2 (D) + r n||/ 


u n + u 71 - 1 




P h (u n ) + Phiu 71 - 1 ) 


A3) 


II 


(3) 


IIl 2 (X>)- 


For the first term we use (33) to get 


1^(3) | < T n M 11^11^2,00(^^2(25)) +T n ^M\\u n - U n II£,4(25) 

= T n M2 \\ u \\w^(I n ,L^(V)) +Tn^M\\ f d t u{-, t) dt\\ 2 L 4 {D) 

• ' In 

< T 3 M 2 ||rt|| W 2,oo(2 n) L2(25)) + TnC M \\u\\wA°°{I n ,L 2 (V))- 


For the term II ( 3 ) we get in the usual manner 

S n 


|III ( 3 )|=r n ||/ 

sn 


u n + u n 1 


-/ 


P h (u n ) + P h (u n ~ l ) 


\ (|25| 

j IIl 2 (X>) RiCjif ||/l ^11 L°°(/ n ,J/ 2 (25))• 


Combining the estimates for I ( 3) and II (3 ) with (41) yields 

S n s n - 


Is^’.si + CVm) 
WEZ + Ef'Wmv) 


< CmRi (II t2 ' u IIvU 2 ' 00 (/„,L 2 (D)) + l|^ 2 ' u ||L 0o (/ n ,H 2 (7?))) 


(3) 

and hence the final estimate for the Qh -term 


l»(<d 3) ,^ + ^- 1 )L»(75 ) ) 


(42) 


< 


r n (ll^h IIl 2 ( 2 )) + ll^ft 1 |Il 2 (D)) + t uCm (ll' r2 ’w|lvU 2 - 00 (i' n ,L 2 (D)) + ll^ w llL“(/„,_ff 2 (D))) • 

Next, we bound the term E%). It holds 

f . . 

i^(d 4) )iiL 2 P) 


iR<L(ei 4) ),^+^- 1 )i p 


|E^ + EM l2 (25) 


(43) 


E( 


u n + 


-u(-,t))dt ||l 2 ( 25 ) < — T 3 ||L( 9 u'ix)|| 2 / oo( 2 n ,Z / 2 ( 25 )). 


Combining the estimates (39)-(43) with the error identity (27) proves the lemma. 


P. 
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Recall that according to Lemma 4.2 and the Poincare-Friedrichs inequality there exist positive 
constants c E and Ce such that 

ce||Vu||| 2(2?) < \\vf E{p) < C e \\S7v\\ 2 l2{V) for all v € (44) 

Lemma 6.5 (Energy-error estimate for Ef). Consider n > 1 and E= U n — Ph{u n ). Let M denote 
the constant in Lemma 5.2. There exists a constant Cm that only depends on M, T>, d, the data 
L, k, (5, the norm-equivalence constants C E and c E and the stability constants C\y i,«>, C E 2 and C\ 
such that for all r n < ( 2Cm)~ 1 it holds 


\Eh\\ E (v) - (1 _ Cmt ) 1 Hs(x>) + ^M\\hdtu\\ 2 L 2^ InH 2(p-)) 


(45) 


+CmTti {\\ t 2 L(dttu)\\ 2 L oo^ In H i^ v ^ + + 11 hu 11 ^ 00 H 2(x>))^j ■ 

Proof. We proceed analogously to the proof of Lemma |6.4| Starting from the energy error identity 
(28) we obtain the following estimates for the various terms. Using (13) we get 

^((LiE^ + E-- 1 )^^)))] 0 

- 1 - S C L 2 \\f^ ^\\e(v) < Cl^CxtJ \\hdtu\\ L 2 ^ InH 2 ^y (46) 

\\ h h+^h II E(V) 

Second, using \\Ph(u n )\\ w i,oo(p ) < C w i,oo (diam(P) + l)\\'Vu n \\ L oo^ < M (cf. (A7)) we get 

|»(<L(£U + ^- 1 ),P La (^ 2 ))))| 


? Ctf\c M PT v (\\]%\\l m + \\El- l \?Em) + Tn Y. ll't(ft(»*)-^*)IU(I>)|K + £;- 1 |U(C) 

44l 

^ C L 2 


k=n— 1 


” 1 ” ^ ^Sobolev 


-pn—1 II2 


" r n||^||vF 1 ’ 3 (X>)(ll^ll J £;(I)) + 11 £(£>), 


(47) 


I Z? n ll^ 1 II 1 11 ^ \ 

\^h\\E(V) + ll£(X>)J 

^AIlF n " 2 

yfcE 

In the last step we also used the following inequality (based on Sobolev embeddings) which holds 
for any v G H l (fD) 

IIV(ku) \\l 2 (v) < ||(V«;)u)|| i 2 ( I j) + \\k(Vv)\\ l2{v) < ||V/cH^s^)Il' t, llz, 6 (x>) + IMIl oo CD)IMI.H' 1 (x>) 

< (C , soboiev||^||rvi.3(x>) + II k IIl°°(:d)) IMI/trd)- 

(o\ 

For the f n ; we use again that |M|l°°(I„x:e>) < M and ||P/t(u ri )||i / °o(x>) < M in combination with 
f M (z) = |^| 2 ^ for |z| < M. This yields 


^((L^ + E-- 1 )^^^))) 

ll^ + K _ 1 lk©) 

< c i2 ^nv 


<C l 2 ^||V^ 3 )|| l2(C) 


(48) 


y/CE 


«(«(■, f) - ( 


PftK) + PftK" 1 ) 


dt ) ||l 2 (x>) 




+/3C i2 ^||V 

VCE 


din 


P h {u n ) + Phiu"- 1 ) 


) dt ] ||l 2 (d)- 


II 


(3) 


2 
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The regularity of k (and the fact that we can hide HrellwMm) and ||«||z,<^>(x>) in Cm) allows us to 
estimate the first term by 1^(3) < CMT n {\\T 2 dttu\\ L oo^i nH np^ + \\hu\\ L oo^ InH 2 (p))). For the second 
term we get 


II f (3) </3C- i2 ^| |vQf fW;t))-f( un+ f V ) ||l 2 P ) 


+r n /3C i 2 ^||V(/( 
\/Fe 


h(3) 

*>n 

u n + u n ~ l 




V) 

sn 


where we can use (34) to obtain 

i^(3) < T^C'Af||u|| w -i,oo(/ Ti ,_ffi(x>)) + ' r nC , M||'ti||ry2,°o(j n ,i : fi(r))) 

and where we can use ||M n ||wi,oo(x>), ||P/ l (w Tl )||^yi,oo(x>) < M to get 

11 ^( 3 ) < CmTti 11 hu | #2 (x))). 


Combining the estimates for I ( 3 ), II ( 3 ), i,( 3 ) and ii ( 3 ) with (48) yields 

sn Sn sn Sn 1—1 


|S({i(E; + Er 1 ),F ia «< 31 ))) 

W+KPTm 


< CmTti (||'T 2 rt|| W 2 ’°°(I n , H l(p)) + \\hu\\L°°(i n ,H 2 (v))) ■ (49) 


For the last term in the error identity (28) we get 


— In r Pt 2 (. l ^)+P t ;^u— 1 ) 


I* [Mez + et 1 ),/, 


- P L 2 (Lu(-,t)) dt}) 


5 ML 


I Jpn i T7' n ~ 1II 
\ h h+ h h II E(V) 


(50) 


n i r„,n— 1 


Lu n + Lu 


- Lu(-,t) dt I || E(V) < T nCL 2 \\L(dttu)\\L°°(I n ,E(T>))- 


Combining estimates (46)-(50) and plugging them into the error identity (28) finishes the proof. □ 
Lemma 6.6 (Full L 2 -error estimate for EC). We use the notation and the assumptions of Lemma 


6.4 Then it holds 


I TP n II z ^ ^ 

\ tj h IIl 2 (X>) — e 


ip: 


°lli 2 (x>) + CMe A ° Mtn \\h 2 dtu\\ 2 L 2^ 0 tritH 2 (p^ 


11 

+CMe. ACAltn (ll r2 F(3tttt)||2oo(/ fe)i 2(X))) + ll r2tt lliV 2 ' 00 (7 fe ,L 2 (I’)) + II^ 2 ' U ||2°°(4,J7 2 (D))) • 


k =1 


Proof. First we note that if a n ,b n ,a n is a sequence of positive real numbers that is related via 
a n+i < (1 + a n )a n + b n then it holds 


&n +1 S: 


< e Er= 0 <*i | Q0 + ^ j. 


(51) 


i =0 


Next we use equation (38) to obtain 

7i n ii2 


E h\\L 2 {T>) ^ (l + «n)|| llz,2(x>) + &n> 


(52) 
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where a n := ( ^ n Tn) and 


b n C M \\h dtu\\ L 2( IritH 2(D)) 

+ CmTu L(dttu) llz / 00 (/ rl) L 2 (X>)) + \\ t2u \\ W 2 ’°°(I n ,L 2 (T>)) + \\h 2u \\ L°°(I n ,i? 2 (X>)) 


Combining (52) with (51) and Cm^u < 1/2 yields 


e < \\E®\\ 2 L 2 (p) + Cm ^2 

fc=i 

n 

+CmE T * (ll r "^'(^ M )lli 00 (4,L 2 (X>)) + ll r2 ' u llv7 2 ’ 00 (I fc ,L 2 (X>)) + \\h~ u \\ 2 L°°(I k ,H 2 (D)) 


k =1 


The inequality Ylk =l — 4 Cm 4 finishes the proof. 


□ 


Lemma 6.7 (Full energy-error estimate for E?). We use the notation and the assumptions of 
Lemma 6.5). It holds 


IK 


n || 2 


7ill£(Z>) ^ 


< e 




\E ' 


'Oil 2 | /O ■ 


iCMtn ll^ u lli 2 ( 0 ,t„,.ff 2 (X>)) 


n 

+CMe 4CMtn ^2 T k (ll r2 -^(^ M )llio°(4,_H' 1 (D)) + \\ T2u \\w 2 ’ 00 (I k ,H 1 (V)) + ll^ M ll|o°(7 fe ,_H' 2 (D))) ■ 


fc=l 


Proof. The proof is analogous to the proof of Lemma 6.6 by combining equation (51) with Lemma 
16.51 □ 

Following the ideas of |38] , we want to show that the solution u ^ of the original discrete problem 
0 is identical to the solution Uf of the auxiliary problem ([22]) implying that the estimates in 
Lemma 6.6 and 6.7 hold equally for u 1 f l . For that purpose, we want to show that if r n is sufficiently 
small it holds WU^W^tv) < M for all n > 0. Then, by the properties of /m, we obtain equality 
of and UJf. To show the desired boundedness we can use again Lemma 5.1, which guarantees 
IKIIl°°(z>) < CMVvhWl^ for all v h G S h . 

Conclusion 6.8. Let assumptions (A1)-(A7) be fulfilled and let h and r n be such that 4(^ niax + 
r^) —» 0 for h,r n —* 0. Then, for all small enough h and r n , the corresponding solution Uf (i.e. the 
solution for /m as specified in Lemma 5.2) fulfills 

\m\ L °o {v) <M. 


Proof. We have LTjf = E'f + P h (u n ). Using (44) and Lemma 5.1 we get 

\\K\\ L ~<p) < \\P h (u n )\\ L ~ {v) +C^^£ h \\EZ\\ E{ vy 


The term Ph(u n ) is uniformly bounded by (A7) and the Poincare-Friedrichs inequality with \\Ph(u 
C' w i,oodiam(T ) )||Vrt n ||ioo(x)). Let us hence consider the second term. Fixing the model problem (and 
assuming (Al)-(A7)), the only variables are h and r n . With this, we can write Lemma 6.7 as: there 
exists a constant C(M), which is independent of h and r n such that 

\\Eh\\ E {v) < C(M)(h max + max 7-2). 

1 <n<N 

Consequently, for each given e > 0, we can pick h(M) and r n (M) small enough so that 
l h {M)\\El\\ E{D) < C(M)4(M)(h max (M) + max r n 2 (M) < e. 

l<n<N 

If we choose h and r n small enough so that CocV^K^KKIKh \\e(V) < ll' u llw 1 .°°(/„,w 1 ’ oo (X>))) then 
we obtain ||E/£||l°°(x>) < M as desired. Qi 


L°°(V) < 
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Observe that Conclusion 6.8 proves Proposition |3.4| We are now prepared to conclude the proof 
of Theorem 13.51 


Proof of Theorem 3.5. We pick h{M) and r n (M) small enough so that the bound in Conclusion 6.8 
holds true. Since ||E7 £||l°o(£>) < M we obtain from the properties of /m (see Lemma 5.2) that Uf 
must be identical to the solution vTf of for every time step n > 1. Hence, we obtain the splitting 

v%-u n = E% + (I\(u n )-u n ), 


be estimated in the usual matter. 
uq £ H 2 (V ) concludes the proof. 


where E£ can be estimated by Lemma 6.6, respectively Lemma 6.7 and where ( Ph{u n ) — 


can 


A Lagrange-interpolation error estimate for the initial value 

M 


7 Numerical experiments 

In this section we investigate the performance of the one-stage Gauss-Legendre implicit Runge- 


Kutta scheme stated in Definition 3.3 and compare it with the approximations obtained with the 


Backward-Euler method (12) to stress the importance of the discrete mass conservation. We consider 
the computational domain V := [—6,6] 2 and the time interval [0, T max ] := [0,100]. We seek a solution 
u £ C°([0,T), Hq(T>)) to the time-dependent Gross-Pitaevskii equation 


idtu = —-A u + V u — 8lC z u + /3\u\ 2 u 


in T>, 


(53) 


where we recall C z = — i (xd y — yd x ). We use the following configuration. We chose /3 = 100, 
f! = 0.8 and the harmonic potential 


P(x) : = 


3x x2 + TyV 2 


with trapping frequencies 7 X = 0.9 and 7 X = 1.1. The initial value uq = u(-, 0) is chosen as the 
L 2 -normalized ground state eigenvector of the Gross-Pitaevskii operator Gq(v) := —^Av + Vq v — 
0.8 C z v + 100|u| 2 u with Vo( x ) = \( x2 + V 2 ) an d corresponding ground state energy Eq = 3.1938 
(cf. [IS]). We computed this ground state using the Discrete Normalized Gradient Flow method 
proposed in m- Starting from this setting, we wish to simulate the dynamics of no in the anisotropic 
trap V, i.e. we solve (53) numerically. The problem is picked in such a way that vortices, i.e. density 
singularities, are formed in the condensate (see Figure [3] and [l]). We define the energy by 


8 (y) := (v,v) E ( V) + 




/ | i2 i ^ \ |4\ 

(k\v\ + — |u| ' 


P\ 

2 


which is a conservative property of equation (53). 


We demonstrate the efficiency of the Crank-Nicolson-type IRK scheme (as stated in Definition 
3.3) by showing that large time steps are allowed, thanks to the mass conservation property. The 


Backward-Euler approach on the other hand (despite being unconditionally stable) does not allow 
large time steps, since this results in a severe loss of mass which lets the corresponding approxima¬ 
tions vanish quickly. In all our computations we use a uniformly refined triangular mesh Th with 
66.049 nodes. That means that the discrete space Sh contains 132.098 degrees of freedom (minus 
the ones from the boundary condition). We use uniform time steps and denote r := r n for simplicity. 

We note that the computational complexity of the IRK scheme ( |11[ ) and the Backward-Euler 
scheme (12) is roughly the same in our implementation. Since both schemes are implicit, they require 


an iterative Newton method in each time step where we observe a comparable number of iterations 
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t=60 t=80 t=100 



Figure 1: The figures show the approximations for the particle density obtained with the Crank- 
Nicolson scheme for large time steps r = 0.1 at times t = 0, t = 20, t = 40, t = 60, t = 80 and 
t = 100. We observe a reduction of the number of vortices. The mass is fully preserved and the 
energy up to a relative error of 0.026%. 

Table 1: The table shows the loss in energy and mass of the Backward-Euler scheme after 100 time 
steps for different time step sizes r, compared wth the corresponding quantities obtained with the 
Crank-Nicolson-type IRK scheme. Recall that the IRK always conserves the mass. 


T 

T 

W u h,BE\\ l 2 (V) 

W U h,IRK\\ l 2 (V) 

£( u h,BE) 

£( u h,IRI<) 

l 

100 

2.6- 10 -17 

1.0 

6.9 • 10 -34 

3.36455 

lO” 1 

10 

0.13275 

1.0 

0.02309 

3.19203 

10 -2 

1 

0.91676 

1.0 

2.52235 

3.19383 

10 -3 

0.1 

0.99905 

1.0 

3.18549 

3.19383 


to reach a given tolerance. In the following, we shall denote Backward-Euler approximations by 
u hBE an( I IRK approximations by u r f j RK . 

Due to the structure of the problem, we could use exact integration when assembling the system 
matrices and load vector for our problem. Furthermore, all linear systems were solved with an 
UMFPACK direct solver. The only reason why we were not computationally exact (up to machine 
precision), was that we prescribed a residual tolerance of order 0(1O~ 8 ) for the Newton-algorithm 
to abort. This inexactness did not have an observable effect on the conservation of mass for the 
IRK in any of the computations. Concerning the energy, a small deviation from the exact energy 
was observable over time for the IRK, however in a negligible range. For instance, for large steps 
r = 1, the energy was still preserved up to an error of 5.3% at T = 100. For slightly smaller 



















24 



Figure 2: The figures shows the approximations for the particle density obtained with the Backward- 
Euler (upper row) and the IRK scheme (lower row) for the time step size r =1 at times t = 20 and 


t = 40. 




Figure 3: The figure depicts numerical approximations for the particle density p(-,T) := \u(-,T)\ 2 
at T = 1. The left approximation is a Backward-Euler approximation computed with time step size 
t = 10“ 2 , whereas the IRK approximation on the right is computed in one time step, i.e. with r = 1. 

time steps with t = 0.1 the conservation of energy already improved to a relative error of below 
0.03% which is insignificant considering that the reference energy (at t = 0) is typically already 
polluted by discretization errors. Using a time step size r = 0.1 we simulated the dynamics of the 
particle density on the time interval [0,100]. The corresponding results are depicted in Figure [I] We 
observe that the condensate with initially seven vortices collapses to a condensate with six vortices 
at T = 100. 

This is in strong contrast to the Backward-Euler scheme that, though unconditionally stable, 
suffers from a major loss of energy and mass. This is clearly shown in Table [l] For time step sizes of 
order r = 1, basically all energy and mass is lost after 100 time steps. In our example the situation 
gradually improves with decreasing time steps sizes, however, to obtain an acceptable loss of mass 
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and energy after 100 time steps, the Backward-Euler method requires time steps sizes of at least 
t = 10 -3 . The significance of the preservation properties is further emphasized in Figure [ 2 ] Here 
we compare Backward-Euler and IRK approximations for large time steps r = 1. We can see that 
even though \v% IRK \ 2 is not particularly accurate, it still preserves the structure of the condensate, 
whereas Iu^beI 2 quickly collapses into a vanishing mass that fully contradicts the correct physical 
behavior. 

Finally, in Figure [3] we compare the IRK approximation after one single step of order t = 1 
with the Backward-Euler approximation at the same time (T = 1) but using 100 time steps with 
size r = 10~ 2 each. Even though the costs for the Backward-Euler scheme are a 100 times higher 
as for the IRK approach to obtain a comparable result, the approximation IvJj, BE \ 2 has clearly not 
yet the quality of \u][ IRK \ 2 . 

In summary we conclude that even though the Backward-Euler scheme seems to be uncondi¬ 
tionally stable, the loss of mass and energy has a tremendous impact on the quality of the obtained 
approximations if the time-step size is not chosen very small. The IRK scheme of Crank-Nicolson- 
type on the other hand does not appear to have such restrictions. 
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